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ABSTRACT 

The inverse computation for the design of phase masks in two layers is investigated and a novel algorithm 
presented. Phase masks are also referred to as phase-only holograms or difiractive optic elements (DOE). A 
two layer system can transform an arbitrary specified light field at an input plane to a desired light field at an 
output plane. The light field includes both intensity and phase. Such a system can be cascaded for higher level 
functionality. There are two computations involved. The first computes a sensitivity matrix symbolically. The 
elements of the matrix hold the variation in each element at the output plane with variation in each element of 
both phase screens. An element of this matrix is provided for reference. The second algorithm iteratively updates 
the phase screen values to bring the output fieldto that desired. On each iteration, the algorithm performs a 
forward computation from input to output. The phase values are updated using the sensitivity matrix and the 
error at the output relative to that desired. A unique solution is not required, only one that provides the required 
transformation from input plane to output plane. 

1. INTRODUCTION 

Two layer diffractive optic element (DOE) or hologram systems, shown in figure 1, are of interest because 
they allow a specified light field in an input plane to be redirected into a desired light field at an output plane. 
The light field has both intensity and phase at every point in the plane transverse to the direction of propagation. 
This allows such units to be cascaded and hence capable of being used as building blocks for more complicated 
structures. Further, changing the propagation model from free space to waveguide allows the same approach for 
designing holographic elements to control light in multi sequence photonic integrated circuits, PICs, or for use 
in interconnecting SOAs for arbitrary high speed logic [4]. Two sequential phase screens may also be used for 
free space interconnection logic [7](Ch. 9). 
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Figure 1: Two layer DOE illustration. 



Section 2 discusses a widely used approach for synthesizing a single phase screen for comparison and then 
describes the two layer synthesis problem. Section 3 provides a summary of the proposed algorithm for computing 
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Figure 2: Synthesis for design of single plane computer generated hologram or DOE. 







Fast Fourier 
transform 
















Phase- 
only 

hologram. 

Intensity 

constant 




Desired 
field 
Intensity 
specified 








1 


diverse fast 
Courier transf orrr 






1 







Figure 3: Gerchberg-Saxton algorithm for single hologram or DOE design. 



phase screens in two layer systems, both forward and inverse computations^ Section 4 and 5 provide details of 
the forward propagation computation and the inverse computation respectively. A conclusion is presented in 
section 6. 

2. CURRENT AND PROPOSED APPROACHES 



Most current approaches to synthesizing holograms or DOEs cannot produce cascadable units. One such 
single DOE algorithms is described for the purpose of illustration in section 2.1. The problem for the two layer 
DOE is described in section 2.2. 

2.1. Phase-only single DOE synthesis approach. 

A number of numerical methods have been developed for synthesizing a single hologram or diflfractive optic 
element (DOE) such that light entering a diffractive optic element is directed into an output plane with a desired 
intensity [11] as shown in figure 2. To minimize loss of light in passing though the hologram, the latter is often 
restricted to phase only. A lens is used to focus light to the output plane placed at the focal point of the lens. 
This allows propagation to be modeled by a Fourier transform. 

A typical widely used algorithm in this case is that by Gerchberg and Saxton [1] as shown in figure 3. A 
plane wave entering the DOE is assumed. The propagation between DOE and output plane is modeled by taking 
a Fourier transform (FT) of the DOE output. At the output plane, the intensity is changed to that desired at 
the output plane. The phase is left at whatever value it has, as there is no requirement on output plane phase. 
An inverse Fourier transform (IFT) is used to transform the signal back to the DOE plane. At this point the 
intensity is set to one as we desire only a phase variation. Repeated iterations of this loop will generate the 



optimum phase DOE to produce the required output intensity at the output plane. As only intensity is specified 
at the output plane, such a system cannot be cascaded. 

A generalization of this system excludes the lens, allowing the lens function to be performed by the hologram, 
and replaces the Fourier transform by a Fesnel diffraction formula. This allows the distance between hologram 
and output plane to be varied within the Fresnel diffraction zone, normally from millimeters to several meters. 

2.2. Two layer DOE synthesis problem 

In this paper we examine the system shown in figure 1. We assume light of specified intensity and phase at 
an input plane that diffracts into the two layer system. Because of the similarity with layered models used for 
the propagation of light through turbulence [6] [5], we use the name phase screen for the phase-only DOE or 
phase-only hologram. In the case of turbulence, the effects of turbulence are localized in phase screens and light 
is propagated between phase screens as though turbulence was absent. In our case, the light also propagates 
through free space between planes. The light propagates from the input plane to the first phase screen. It then 
passes through the first phase screen and propagates to the second phase screen. It passes through the second 
phase screen and propagates to the output plane. The phase screen phase values are computed to produce a 
specified intensity and phase at the system output plane, given a specified intensity and phase field at the system 
input plane. 

In figure 1, we also show a way to computationally cope with the diffraction spreading into larger transverse 
areas. The second phase screen is padded around with zeroes on all sides so as to transform to a larger area 
output. 



3. SUMMARY OF ALGORITHM FOR COMPUTING PHASE SCREENS IN TWO LAYER 

SYSTEM 

3*1* Forward computation 

All operators and fields are assume to be functions of orthogonal directions x and y transverse to the direction 
of propagation z. Light fields are assumed to have both intensity and phase defined at every x — y element in 
the field. While transverse fields are continuous, we will discretize them for the purpose of computation, with 
sufficient resolution to maintain the highest spatial frequencies desired. The following is a summary, the individual 
items are elaborated subsequently. 

Propagation of the light field from one plane to another is described by the operator D. The spatially 
discretized form of D will be a matrix. Propagation through a phase screen is accomplished by modifying each 
element of the incoming light by adjustment of phase <f>. If the incoming 2-D plane of light is rearranged as a 
vector, element by element multiplication for the phase screen can be represented by multiplication by a diagonal 
matrix H with elements ff n , m = exp{j<f> nim }. As we plan to determine the optimum H for each phase screen, 
we use small random values for the first iteration or an educated guess and then iterate to improve the values 
so as to meet the required output field. Further, we use a subscript to indicate which phase screen or region is 
in consideration starting with 1 at the left in figure 1. Therefore, for the ith iteration and the jth phase screen 
in the layered model we write Hj(i). The forward computation from input plane to output plane at the ith 
iteration can then be represented by: 

F(0 = -D 3 H 2 (0^2ffi(0^i (1) 
Given an input field X y the field in the output plane can be written: 

y(0 = A>#2(WM0£i* (2) 

For convenience in subsequent steps the two phase screens are included as partitions into a single matrix: 



*(»■) = [Mfc(Q] 



(3) 



As an example, used for computing the sensitivity element shown in the appendix, for two 1-D four long phase 
screens the phase values can be written: 

$(i) = [pll pl2 pl3 P U\p21 p22 p23 p2A ) (4) 

Note that, figure 1, as the output is twice the size of the phase screens, there are the same number of unknown 
phase terms in the equations as there are terms in the output plane. 

3.2. Inverse computation 

A square sensitivity matrix S(i) is computed using: 



The output field at this point will differ from the desired output field Y4. Therefore we compute difference 
vector Ay (i) for every element in the 2-D plane. 

AK(t) = r(i)-r, (6) 

The phase screen values for the ith iteration are substituted into the sensitivity matrix and the incremental 
update for the phase screens A$ is obtained by solving the equation: 

5(i)A*(i) = Ar(i) (7) 

For realistic size problems the ill conditioning of S(i) must be taken into account in finding an acceptable solution. 
The phase screens are now updated for the next iteration: 

*(i + l) = *(0 + A*(i) (8) 

We now return to equation (2) with i + 1 in equation (6) replaced by i to perform the next iteration. We stop 
iterating when the normalized changes are sufficiently small, for example: 

lA^fOI 2 

, < desired accuracy (9) 

We also stop iterating if the error no longer decreases noticeably from one iteration to the next. In this case we 
may need to change the phase screen size or sampling to come closer to the desired accuracy. 

4. DETAILS OF FORWARD COMPUTATION FOR COMPUTING PHASE SCREENS IN A 

TWO LAYER SYSTEM 

4.1. Details of propagation in free space between planes 

The propagation distance between planes , z, is chosen constant for conveniences and is assumed much 
greater than the maximum dimensions transversely in x and y. In this case the Helmholtz wave equations can 
be approximated by the paraxial wave equation: 

Propagation from one plane to another across the layer, across z, can be accomplished by using the Huygens- 
Fresnel principle [2]. This Frensel diffraction can be written as a convolution integral for homogeneous media: 



V(x, y) = / U(i,f})h(x - C,» - eta)dtdT) 



(11) 



where the impulse response is: 



The transfer function is the Fourier transform of the impulse response [2]: 



(12) 



(13) 



Note that convolution in the space domain becomes multiplication in the Fourier transform domain. Therefore 
propagation of a field U through a distance Az is achieved by taking its Fourier transform, multiplying by the 
transform function in equation (13) and then taking the inverse Fourier transform: 



V{x, y, z + Ax) = ?- 1 {TV{x, y, z)H{k Xi ky)} 



(14) 



4.2. Details of propagation through a phase screen 

hi propagating through a phase screen, each element of the field in the transverse x-y plane is multiplied 
by exp{<j>(i)}. If the 2-D transverse field is •written as a vector, the operation can be written as multiplication 
by a diagonal matrix. For example a 2 x 2 phase screen <j>(n,m) givesa matrix: 



H = 



expifarf)} 0 0 0 

0 expifod)} 0 0 

0 0 exp{<h,i{x)} 0 

0 0 0 eip{02,2(i)} 



(15) 



where the diagonal elements are the transverse values exp{<f> n , m (i)} at the transverse x-y coordinate for the 
phase screen at the ith iteration. 

5. DETAILS OF COMPUTATION OF THE SENSITIVITY MATRIX AND INVERSE 

COMPUTATION 

5.1, Symbolic propagation through a phase screen with phase values as variables 

In order to compute a sensitivity matrix we perform a symbolic discrete Fourier transform for propagation 
through free space instead of a Past Fourier transform as used in section 2.1. For example, a 1-D discrete Fourier 
transform from u to U for vectors of length N = 4 can be written as a matrix-vector multiplication: 
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where dx = A/10 is the interval between samples in the space domain, A is the wavelength of the light and 

W = exp{2**j/N) (17) 
Similarly, the 1-D discrete inverse Fourier transform from V to u can be written as a matrix-vector multiplication: 
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(18) 



where dfx is the interval between samples in the spatial frequency domain. Note that dx * dfx = 1/N. Matrix 
multiplication of equation (16) by equation (18) gives the identity matrix, consistent with forward and inverse 



transformations canceling ; each other out. Note that the Fourier transform and inverse Fourier transform per 
formed for D 3 h^ twice the dimension of that for D, and D 2 as explained previously. For 2-D phaTscrSn £ 
need to use 2-D Fourier transforms. eens we 

5.2. Symbolic computation of sensitivity matrix with phase values as variables 

^ Performing the forward confutation symbolically with variables for the phase values allows us to make the 
drferentiation m equation (5). In the appendix we show the first element out of 64 for an 8 x 8 sensitivity matrix 
obtained for four long phase screens. As can be seen the sensitivity matrix is a lengthy equation for any practS 
system. However, as shown m an earlier reference [9] [10], the results can be reformulated to a more impact 

5.3. Inverse computation to estimated phase values in screens 

Eq^tion (7) is a nonlinear equation for the phase value corrections because the sensitivity matrix depends 
on the iteration and current phase. It is solved by linearization in the iteration method proposed. Fo^W 
reahstic matrices, Gauss elimination and similar methods take too long. The most practical technique for solviS 

&b£ ^."^WftaS * ™>nna% » * the conjugate gradient method 
Solution time depends on the clustenng of eigenvalues in the sensitivity matrix after the current phase values 
have been mserted Further t is normal to condition the matrix, first by normalizing rows and hZ^Z 
tte normalization at a later tune and second by preconditioning the matrix to improve the eigenstructure The 
preconchtaon must also be backed out later. The conjugate gradient method will provide the minimum en^ 
SmS? ^."ormaUy be perfectly adequate to achieve the goal of correctly directing the intensity^ 
phase of the hght ink .the output plane. It is not required that a unique solution be found, rather any that 
adequately performs the necessary function will suffice. We are currently validating the approach for design of 
two layer phase screen systems and have not yet adequately verfied the result to warrant inclusion. 

6. CONCLUSION 

pwr;^^ ° f deSigDing P ^ C SCTeenS ' *° ^ 0wn 38 P^e-only holograms or dim-active optical 
tZTv ?r,^° + B f re ! ns 1 are used * S «W«™ to change a specified hght field at an input plane to a 
desn-ed hght field at an output plane. In this case the intensity and phase at the input plane and the desired 
mtenaty and phase a the output plane are fully specified. The method allows the output plane to be 
than the mput plane for diffiaction spreading. Such units can be cascaded to create more complex structures 
unlike the case for a single hologram or DOE. H 

r^LT^t^f^ 311(1 5 DUmeriC ^P^tion. In the symbolic computation, the forward 

compu ation of the field from input to output is performed using discrete Fourier transforms and leaving phase 

t^Z l + K I SCT T "J"*** A S6nSitivity matrbc 13 COm P uted b y differentiation of the output with 
SwSL T P t? demeDtS repreSent the «*■»»■ * out Put elements with respect to change 

m phase value. The symbohc first element of the matrix is recorded for reference. The numerical computation 
performs the following . operations rteratively until the output field adequately represents the desired field. With 

S^TTv^S P ^ aS - e ' ft ^ rWard com P u t ati on * computed to obtain an output field. The difference 

betweer .this field and the desired field is used with the sensitivity matrix to compute an update to the phase 
values. The algorithm presented is very computationally demanding and is practical only because of the increased 
T^JST* T r ^ lmpr ° Ved devel °P ment s in handling large matrix ill-conditioned equations. 

SlSkT"!.?!^ reqmred ' 0nly ° De that ^P^ 5 the goal of converting hght in the input 
plane to the desired field at the output plane. 

In addition to verifying the algorithm for specific cases, an important extension is possible. Replacing the free 
cTcSt^pS) ^ * Pr ° Pagati0n ^ the d6Sign ° f cascad able logic as in a photonic integrated 
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Appendix A 



The first element of an 8 x 8 sensitivity matrix, for two four long 1-D phase screens, is shown for illustrative 
purposes and to provide a reference for anyone else attempting this approach. pll,pl2,pl3, jpl4, j?21,p22,p23,p24 
represent the phase elements for the two phase screens. 

> elementli:=evalf (sens2[l,l] ,3) ; element 11 :« 

> (- . 418e6- . 721e6*D* (( . 883e-8- • 612e-7*I)*p21- ( .265e-7+ .306e-7*I) *p22+ ( . 

> 265e-7+ . 306e-7*I) *p23- ( . 265e-7+ . 306e-7*I ) *p24) + ( . 833e6-753 . *I) * ( ( - . 6 12 

> e : 7- • 883e-8*I) *p21+ (-.. 290e-8+ . 404e-7*I) *p22- ( . 265e-7+ . 306e-7*I) *p23+ ( ; 

> 404e-7+ . 290e-8*I) *p24) - ( . 417e6+ .722e6*I) * ( (- • 883e-8+ . 612e-7*I) *p21+ ( . 3 

> 06e-7-.265e-7*I)*p22+(.266e-7+.306e-7*I)^23+(-.306e-7+.265e-7*I)*p24) 

> - ( . 417e6+ . 722e6*I ) * ( ( . 612e-7+ . 883e-8*I) *p21- ( . 404e-7+ . 290e-8*I) *p22- ( . 

> 265e-7+ . 306e-7*I) *p23+ ( . 290e-8- .404e-7*I) *p24) + ( . 833e6+0 . *I) * ( ( . 883e-8 

> ~.612e-7*I)^21+(.265e-7+.306e-7*I)*^ 

> -7+.306e-7*I)*p24)-(.417e6+.722e6*I)*(^^ 

> 8-.404e-7*I)*p22-(.265e-7+.306e-7*I)*p23-(.404e-7+.290e-8*I)*p24)-(.41 

> 7e6+,722e6*I)*((-.883e-8+.612e-7*I)*p21+(«-.306e-7+.265e-7*I)*p22+(.265 

> e-7+ . 306e-7*I) *p23+ ( . 306e-7- . 265e-7*I) *p24) + ( . 833e6-753 . *I) * « . 612e-7+ 

> .883e-8*I)*p21+(.404e-7+.290e-8*I)^22-(.265e-7+.306e~7*I)*p23+(-.29Oe 

> -8+.404e-7*I)*p24) 
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